library(rgl)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
このTDAパッケージはパーシステントホモロジーのライブラリとしてあちこちに見られるもの(GUDHI,Dionysus,PHAT)をベースにしたRラッパーらしいので、これを使うのは多分正解。
このパッケージの Vigneteをなぞろう。
library(TDA)
二次元平面に円があるとして、そこから取られた乱点を作る。 これが観測データとして、パーシステントホモロジーを使って、「形状」などについて解析を加えていく。
X <- circleUnif(400)
plot(X)
次に、この円状オブジェクトの置かれている空間を評価する為に、長方形領域のグリッド点を発生させる。
Xlim <- c(-1.6,1.6)
Ylim <- c(-1.7,1.7)
by <- 0.065
Xseq <- seq(from=Xlim[1],to=Xlim[2],by=by)
Yseq <- seq(from=Ylim[1],to=Ylim[2],by=by)
Grid <- expand.grid(Xseq,Yseq)
plot(Grid,pch=20,cex=0.1)
グリッド点のそれぞれについて、円周状の点集合からの「距離」を計算する。 このdistFct()関数は、Xの点の中でもっとも短いユークリッド距離の二乗(二乗ノルム)の最小値を返すらしい。
distance <- distFct(X = X, Grid=Grid)
DD <- distance
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
別の方法で、グリッド点と円周状の点集合からの「距離」とを計算する。 DTMはm0をパラメタとしてスムージングしながら「距離」を定める仕様になっている。 DTMの定義は、ビグネット https://cran.r-project.org/web/packages/TDA/vignettes/article.pdf を参照。
m0 <- 0.1
DTM <- dtm(X=X,Grid=Grid,m0=m0)
DD <- DTM
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
m0値を変えてみよう。
m0 <- 0.99
DTM <- dtm(X=X,Grid=Grid,m0=m0)
DD <- DTM
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
ここまではグリッド点と点集合Xとの距離を計算していた。
点集合から、その背景にある「分布」を推定してみる。
k-nearest neighbor 法でグリッド各点の確率密度を推定している。 kの値が小さければ細かく、大きければ粗く推定される。
k <- 60
kNN <- knnDE(X=X,Grid=Grid,k=k)
DD <- kNN
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
k <- 100
kNN <- knnDE(X=X,Grid=Grid,k=k)
DD <- kNN
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
ガウシアンカーネルでの密度推定もできる。 hはその粗さを決めるパラメタ
h <- 0.3
KDE <- kde( X = X, Grid=Grid,h=h)
DD <- KDE
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
小刻みにしてみる。
h <- 0.1
KDE <- kde( X = X, Grid=Grid,h=h)
DD <- KDE
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = "lightblue", border = NA, d = 1, scale = FALSE,expand = 0.5, shade = 0.9,ticktype="detailed")
image(matrix(DD,length(Xseq),length(Yseq)))
グリッド点に距離を定めることと、密度推定することとに関係があることがわかったので、ガウシアンカーネルを用いて、密度推定する代わりに、「距離」を算出することにする。
h <- 0.3
Kdist <- kernelDist(X=X,Grid=Grid,h=h)
DD <- Kdist
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
h <- 0.1
Kdist <- kernelDist(X=X,Grid=Grid,h=h)
DD <- Kdist
persp(Xseq, Yseq,matrix(DD, ncol = length(Yseq), nrow = length(Xseq)), xlab = "",ylab = "", zlab = "", theta = -20, phi = 35, ltheta = 50,col = 2, border = NA, main = "", d = 0.5, scale = FALSE,expand = 3, shade = 0.9)
image(matrix(DD,length(Xseq),length(Yseq)))
ブートストラップ法を用いて、グリッド点に与える値(距離・密度)の区間推定をする関数も実装されている。
band.dtm <- bootstrapBand(X = X, FUN = dtm, Grid = Grid, B = 100, parallel = FALSE, alpha = 0.1, m0 = 0.1)
band.knnDE <- bootstrapBand(X = X, FUN = knnDE, Grid = Grid, B = 100, parallel = FALSE, alpha = 0.1, k=60)
band.kde <- bootstrapBand(X = X, FUN = kde, Grid = Grid, B = 100, parallel = FALSE, alpha = 0.1, h=0.3)
繋がったり、穴が潰れたりの様子を計算して、その結果をプロットする。
DiagGrid.kde <- gridDiag(X = X, FUN = kde, h = 0.3, lim = cbind(Xlim, Ylim), by = by,sublevel = FALSE, library = "Dionysus", location = TRUE,printProgress = FALSE)
plot(DiagGrid.kde[["diagram"]], band = 2 * band.kde[["width"]])
DiagGrid.dtm <- gridDiag(X = X, FUN = dtm, m0 = 0.1, lim = cbind(Xlim, Ylim), by = by,sublevel = FALSE, library = "Dionysus", location = TRUE,printProgress = FALSE)
plot(DiagGrid.dtm[["diagram"]], band = 2 * band.dtm[["width"]])
結果の表示を変える。
par(mfrow = c(1, 2), mai = c(0.8, 0.8, 0.3, 0.1))
plot(DiagGrid.kde[["diagram"]], rotated = TRUE, band = band.kde[["width"]],main = "Rotated Diagram")
plot(DiagGrid.kde[["diagram"]], barcode = TRUE, main = "Barcode")
Circle1 <- circleUnif(60)
Circle2 <- circleUnif(60, r = 2) + 3
Circles <- rbind(Circle1, Circle2)
plot(Circles)
maxdimension=1は「形」が1次元多様体であることを指定(か?)
maxscale <- 5 # limit of the filtration
maxdimension <- 1 # components and loops
DiagRips <- ripsDiag(X = Circles, maxdimension, maxscale,library = c("GUDHI", "Dionysus"), location = TRUE, printProgress = FALSE)
plot(DiagRips[["diagram"]])
形を単純に戻してパーシステントホモロジー解析的な(Filtrationによる)、「形」の答えを出して描く。計算時間を短くする為に、少数の点にして、まずやってみる。
X <- circleUnif(n = 30)
# persistence diagram of alpha complex
DiagAlphaCmplx <- alphaComplexDiag(X = X, library = c("GUDHI", "Dionysus"), location = TRUE,printProgress = TRUE)
## # Generated complex of size: 115
##
## 0% 10 20 30 40 50 60 70 80 90 100%
## |----|----|----|----|----|----|----|----|----|----|
## ***************************************************
## # Persistence timer: Elapsed time [ 0.000165 ] seconds
# plot
par(mfrow = c(1, 2))
plot(DiagAlphaCmplx[["diagram"]], main = "Alpha complex persistence diagram")
one <- which(DiagAlphaCmplx[["diagram"]][, 1] == 1)
one <- one[which.max(DiagAlphaCmplx[["diagram"]][one, 3] - DiagAlphaCmplx[["diagram"]][one, 2])]
plot(X, col = 1, main = "Representative loop")
for (i in seq(along = one)) {
for (j in seq_len(dim(DiagAlphaCmplx[["cycleLocation"]][[one[i]]])[1])) {
lines(DiagAlphaCmplx[["cycleLocation"]][[one[i]]][j, , ], pch = 19,cex = 1, col = i + 1)
}
}
par(mfrow = c(1, 1))
" user-defined filtration“として” simplicial complex and the filtration values on each simplex“。
上で用いた関数のカウンターパートは以下の通り。
ripsDiag —> ripsFiltration
alphaComplexDiag —> alphaComplexFiltration
alphaShapeDiag —> alphaShapeFiltration
円柱状の点集合を作る。
n <- 60
X <- cbind(circleUnif(n = n), runif(n = n, min = -0.1, max = 0.1))
#X <- X + rnorm(length(X),0,0.2)
plot3d(X)
DiagAlphaShape <- alphaShapeDiag(X = X, maxdimension = 2, library = c("GUDHI", "Dionysus"), location = TRUE,
printProgress = TRUE)
## # Generated complex of size: 1463
##
## 0% 10 20 30 40 50 60 70 80 90 100%
## |----|----|----|----|----|----|----|----|----|----|
## ***************************************************
## # Persistence timer: Elapsed time [ 0.006325 ] seconds
You must enable Javascript to view this page properly.
ripsFiltration()関数により、単体的複体が得られる。
maxscale <- 0.4 # limit of the filtration
maxdimension <- 1 # components and loops. 1->線とループ、0->線のみ
FltRips <- ripsFiltration(X = X, maxdimension = maxdimension,maxscale = maxscale, dist = "euclidean", library = "GUDHI",printProgress = TRUE)
## # Generated complex of size: 698
par(mfrow = c(1, 2))
plot(DiagAlphaShape[["diagram"]])
plot(X[, 1:2], col = 2, main = "Representative loop of alpha shape filtration")
one <- which(DiagAlphaShape[["diagram"]][, 1] == 1)
one <- one[which.max(DiagAlphaShape[["diagram"]][one, 3] - DiagAlphaShape[["diagram"]][one, 2])]
for (i in seq(along = one)) {
for (j in seq_len(dim(DiagAlphaShape[["cycleLocation"]][[one[i]]])[1])) {
lines(DiagAlphaShape[["cycleLocation"]][[one[i]]][j, , 1:2], pch = 19,cex = 1, col = i)
}
}
par(mfrow = c(1, 1))
“cmplx”は単体的複体を使うというオプション。
m0 <- 0.1
dtmValues <- dtm(X = X, Grid = X, m0 = m0)
FltFun <- funFiltration(FUNvalues = dtmValues, cmplx = FltRips[["cmplx"]])
DiagFltFun <- filtrationDiag(filtration = FltFun, maxdimension = maxdimension,library = "Dionysus", location = TRUE, printProgress = TRUE)
##
## 0% 10 20 30 40 50 60 70 80 90 100%
## |----|----|----|----|----|----|----|----|----|----|
## ***************************************************
## # Persistence timer: Elapsed time [ 0.001099 ] seconds
par(mfrow = c(1, 2), mai=c(0.8, 0.8, 0.3, 0.3))
plot(X, pch = 16, xlab = "",ylab = "")
plot(DiagFltFun[["diagram"]], diagLim = c(0, 1))
plot(X)
#X <- matrix(rnorm(length(X)),ncol=2)*0.5
for(i in 1:length(FltRips[[1]])){
tmp <- FltRips[[1]][[i]]
if(length(tmp)==1){
}else if(length(tmp)==2){
segments(X[tmp[1],1],X[tmp[1],2],X[tmp[2],1],X[tmp[2],2])
}else{
tmp2 <- c(tmp,tmp[1])
#print(tmp2)
for(j in 1:length(tmp)){
segments(X[tmp2[j],1],X[tmp2[j],2],X[tmp2[j+1],1],X[tmp2[j+1],2])
}
}
}
plot3d(X)
#X <- matrix(rnorm(length(X)),ncol=2)*0.5
for(i in 1:length(FltRips[[1]])){
tmp <- FltRips[[1]][[i]]
if(length(tmp)==1){
}else if(length(tmp)==2){
segments3d(t(matrix(c(X[tmp[1],1],X[tmp[1],2],X[tmp[1],3],X[tmp[2],1],X[tmp[2],2],X[tmp[2],3]),ncol=2)))
}else{
tmp2 <- c(tmp,tmp[1])
#print(tmp2)
for(j in 1:length(tmp)){
segments3d(t(matrix(c(X[tmp2[j],1],X[tmp2[j],2],X[tmp2[j],3],X[tmp2[j+1],1],X[tmp2[j+1],2],X[tmp2[j+1],3]),ncol=2)))
}
}
}
You must enable Javascript to view this page properly.
2方法ある。